<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<title>Kriging Library: sampling.f90 Source File</title>

<link href="tabs.css" rel="stylesheet" type="text/css"/>
<link href="doxygen.css" rel="stylesheet" type="text/css" />

<link href="search/search.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="search/search.js"></script>
<script type="text/javascript">
  $(document).ready(function() { searchBox.OnSelectItem(0); });
</script>
<script src="../mathjax/MathJax.js">
  MathJax.Hub.Config({
    extensions: ["tex2jax.js"],
    jax: ["input/TeX","output/HTML-CSS"],
});
</script>
</head>
<body>
<div id="top"><!-- do not remove this div! -->


<div id="titlearea">
<table cellspacing="0" cellpadding="0">
 <tbody>
 <tr style="height: 56px;">
  
  
  <td style="padding-left: 0.5em;">
   <div id="projectname">Kriging Library
   
   </div>
   
  </td>
  
  
  
 </tr>
 </tbody>
</table>
</div>

<!-- Generated by Doxygen 1.7.5.1 -->
<script type="text/javascript">
var searchBox = new SearchBox("searchBox", "search",false,'Search');
</script>
  <div id="navrow1" class="tabs">
    <ul class="tablist">
      <li><a href="index.html"><span>Main&#160;Page</span></a></li>
      <li><a href="annotated.html"><span>Data&#160;Types&#160;List</span></a></li>
      <li class="current"><a href="files.html"><span>Files</span></a></li>
      <li>
        <div id="MSearchBox" class="MSearchBoxInactive">
        <span class="left">
          <img id="MSearchSelect" src="search/mag_sel.png"
               onmouseover="return searchBox.OnSearchSelectShow()"
               onmouseout="return searchBox.OnSearchSelectHide()"
               alt=""/>
          <input type="text" id="MSearchField" value="Search" accesskey="S"
               onfocus="searchBox.OnSearchFieldFocus(true)" 
               onblur="searchBox.OnSearchFieldFocus(false)" 
               onkeyup="searchBox.OnSearchFieldChange(event)"/>
          </span><span class="right">
            <a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.png" alt=""/></a>
          </span>
        </div>
      </li>
    </ul>
  </div>
  <div id="navrow2" class="tabs2">
    <ul class="tablist">
      <li><a href="files.html"><span>File&#160;List</span></a></li>
      <li><a href="globals.html"><span>File&#160;Members</span></a></li>
    </ul>
  </div>
<div class="header">
  <div class="headertitle">
<div class="title">sampling.f90</div>  </div>
</div>
<div class="contents">
<a href="sampling_8f90.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="comment">!Copyright (C) 2012 Brian A. Lockwood</span>
<a name="l00002"></a>00002 <span class="comment">!</span>
<a name="l00003"></a>00003 <span class="comment">!This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.</span>
<a name="l00004"></a>00004 <span class="comment">!</span>
<a name="l00005"></a>00005 <span class="comment">!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.</span>
<a name="l00006"></a>00006 <span class="comment">!</span>
<a name="l00007"></a>00007 <span class="comment">!You should have received a copy of the GNU General Public License along with</span>
<a name="l00008"></a>00008 <span class="comment">!this program. If not, see &lt;http://www.gnu.org/licenses/&gt;.</span>
<a name="l00009"></a>00009 
<a name="l00010"></a><a class="code" href="classsampling.html">00010</a> <span class="keyword">module</span> <a class="code" href="classsampling.html">sampling</a>
<a name="l00011"></a>00011 <span class="keyword">contains</span>
<a name="l00012"></a><a class="code" href="classsampling.html#a8aa4044081a9abe90d0b88583d2a9908">00012</a>   <span class="keyword">function </span>dinvnorm(p)
<a name="l00013"></a>00013 
<a name="l00014"></a>00014     <span class="keyword">implicit none</span>
<a name="l00015"></a>00015 
<a name="l00016"></a>00016       <span class="keywordtype">real*8</span> :: dinvnorm,p,p_low,p_high
<a name="l00017"></a>00017       <span class="keywordtype">real*8</span> :: a1,a2,a3,a4,a5,a6
<a name="l00018"></a>00018       <span class="keywordtype">real*8</span> :: b1,b2,b3,b4,b5
<a name="l00019"></a>00019       <span class="keywordtype">real*8</span> :: c1,c2,c3,c4,c5,c6
<a name="l00020"></a>00020       <span class="keywordtype">real*8</span> :: d1,d2,d3,d4
<a name="l00021"></a>00021       <span class="keywordtype">real*8</span> :: z,q,r
<a name="l00022"></a>00022       a1=-39.6968302866538
<a name="l00023"></a>00023       a2=220.946098424521
<a name="l00024"></a>00024       a3=-275.928510446969
<a name="l00025"></a>00025       a4=138.357751867269
<a name="l00026"></a>00026       a5=-30.6647980661472
<a name="l00027"></a>00027       a6=2.50662827745924
<a name="l00028"></a>00028       b1=-54.4760987982241
<a name="l00029"></a>00029       b2=161.585836858041
<a name="l00030"></a>00030       b3=-155.698979859887
<a name="l00031"></a>00031       b4=66.8013118877197
<a name="l00032"></a>00032       b5=-13.2806815528857
<a name="l00033"></a>00033       c1=-0.00778489400243029
<a name="l00034"></a>00034       c2=-0.322396458041136
<a name="l00035"></a>00035       c3=-2.40075827716184
<a name="l00036"></a>00036       c4=-2.54973253934373
<a name="l00037"></a>00037       c5=4.37466414146497
<a name="l00038"></a>00038       c6=2.93816398269878
<a name="l00039"></a>00039       d1=0.00778469570904146
<a name="l00040"></a>00040       d2=0.32246712907004
<a name="l00041"></a>00041       d3=2.445134137143
<a name="l00042"></a>00042       d4=3.75440866190742
<a name="l00043"></a>00043       p_low=0.02425
<a name="l00044"></a>00044       p_high=1-p_low
<a name="l00045"></a>00045       <span class="keyword">if</span>(p.lt.p_low) goto 201
<a name="l00046"></a>00046       <span class="keyword">if</span>(p.ge.p_low) goto 301
<a name="l00047"></a>00047 201   q=dsqrt(-2*dlog(p))
<a name="l00048"></a>00048       z=(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6)/((((d1*q+d2)*q+d3)*q+d4)*q+1)
<a name="l00049"></a>00049       goto 204
<a name="l00050"></a>00050 301   <span class="keyword">if</span>((p.ge.p_low).and.(p.le.p_high)) goto 202
<a name="l00051"></a>00051       <span class="keyword">if</span>(p.gt.p_high) goto 302
<a name="l00052"></a>00052 202   q=p-0.5
<a name="l00053"></a>00053       r=q*q
<a name="l00054"></a>00054       z=(((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q/(((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1)
<a name="l00055"></a>00055       goto 204
<a name="l00056"></a>00056 302   <span class="keyword">if</span>((p.gt.p_high).and.(p.lt.1)) goto 203
<a name="l00057"></a>00057 203   q=dsqrt(-2*dlog(1-p))
<a name="l00058"></a>00058       z=-(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6)/((((d1*q+d2)*q+d3)*q+d4)*q+1)
<a name="l00059"></a>00059 204   dinvnorm=z
<a name="l00060"></a>00060 
<a name="l00061"></a>00061       return
<a name="l00062"></a>00062 
<a name="l00063"></a>00063 <span class="keyword">    end function dinvnorm</span>
<a name="l00064"></a>00064 
<a name="l00065"></a>00065 
<a name="l00066"></a><a class="code" href="classsampling.html#a9e8297bff3aca9c1a4b2fa9082e49b63">00066</a>         <span class="keyword">subroutine </span>Latin_Hypercube(ndim,nsam,bound,sample)
<a name="l00067"></a>00067         <span class="keyword">implicit none</span>
<a name="l00068"></a>00068         <span class="keywordtype">integer</span>, <span class="keywordtype">intent(in)</span> :: ndim,nsam
<a name="l00069"></a>00069         <span class="keywordtype">double precision</span>, <span class="keywordtype">dimension(2,ndim)</span>, <span class="keywordtype">intent(in)</span> :: bound
<a name="l00070"></a>00070         <span class="keywordtype">double precision</span>, <span class="keywordtype">dimension(ndim,nsam)</span>, <span class="keywordtype">intent(out)</span> :: sample
<a name="l00071"></a>00071 
<a name="l00072"></a>00072         <span class="keywordtype">integer</span> :: i,j,k,l,key
<a name="l00073"></a>00073         <span class="keywordtype">double precision</span> :: ran,dv01
<a name="l00074"></a>00074         <span class="keywordtype">double precision</span>, <span class="keywordtype">dimension(0:nsam)</span> :: split
<a name="l00075"></a>00075         <span class="keywordtype">integer</span>, <span class="keywordtype">dimension(ndim,nsam)</span> :: iflag_split
<a name="l00076"></a>00076 
<a name="l00077"></a>00077         <span class="keyword">if</span>(nsam.eq.0)return
<a name="l00078"></a>00078         <span class="keyword">if</span>(nsam.lt.0.or.ndim.le.0)stop<span class="stringliteral">&#39;LHS dimension/nsample&#39;</span>
<a name="l00079"></a>00079         sample = 0.d0
<a name="l00080"></a>00080 
<a name="l00081"></a>00081 <span class="comment">! Set Cartesian</span>
<a name="l00082"></a>00082         <span class="keyword">do</span> 100 i=0,nsam
<a name="l00083"></a>00083            split(i) = dble(i)/dble(nsam)
<a name="l00084"></a>00084 100     continue
<a name="l00085"></a>00085 
<a name="l00086"></a>00086 <span class="comment">! Arrangement</span>
<a name="l00087"></a>00087         <span class="keyword">do</span> 200 i=1,nsam
<a name="l00088"></a>00088            <span class="keyword">do</span> 210 j=1,ndim
<a name="l00089"></a>00089 250           continue
<a name="l00090"></a>00090               call random_number(ran)
<a name="l00091"></a>00091               <span class="keyword">do</span> 220 k=1,nsam
<a name="l00092"></a>00092                  <span class="keyword">if</span>(ran.ge.split(k-1).and.ran.le.split(k))<span class="keyword">then</span>
<a name="l00093"></a>00093                    <span class="keyword">do</span> 230 l=1,i-1
<a name="l00094"></a>00094                       <span class="keyword">if</span>(iflag_split(j,l).eq.k)go to 250 <span class="comment">! already used, try again</span>
<a name="l00095"></a>00095 230                continue    <span class="comment">! previous sample loop</span>
<a name="l00096"></a>00096                    iflag_split(j,i) = k
<a name="l00097"></a>00097                    go to 210
<a name="l00098"></a>00098                  <span class="keyword">end if</span>
<a name="l00099"></a>00099 220           continue         <span class="comment">! split loop</span>
<a name="l00100"></a>00100               stop<span class="stringliteral">&#39;out of order at split loop&#39;</span>
<a name="l00101"></a>00101 210        continue            <span class="comment">! dimension loop</span>
<a name="l00102"></a>00102 200     continue               <span class="comment">! sample loop</span>
<a name="l00103"></a>00103         
<a name="l00104"></a>00104 <span class="comment">! Check</span>
<a name="l00105"></a>00105         <span class="keyword">do</span> 300 j=1,ndim
<a name="l00106"></a>00106           <span class="keyword">do</span> 310 i=1,nsam
<a name="l00107"></a>00107              <span class="keyword">if</span>(iflag_split(j,i).le.0.or.iflag_split(j,i).gt.nsam) &amp;
<a name="l00108"></a>00108              stop<span class="stringliteral">&#39;iflag_split&#39;</span>
<a name="l00109"></a>00109              <span class="keyword">do</span> 320 k=i+1,nsam
<a name="l00110"></a>00110                 <span class="keyword">if</span>(iflag_split(j,i).eq.iflag_split(j,k)) &amp;
<a name="l00111"></a>00111                 stop<span class="stringliteral">&#39;overlap of iflag_split&#39;</span>
<a name="l00112"></a>00112 320          continue
<a name="l00113"></a>00113 310       continue
<a name="l00114"></a>00114 300     continue
<a name="l00115"></a>00115 
<a name="l00116"></a>00116 <span class="comment">! Generation</span>
<a name="l00117"></a>00117         <span class="keyword">do</span> 400 i=1,nsam
<a name="l00118"></a>00118           <span class="keyword">do</span> 410 j=1,ndim
<a name="l00119"></a>00119              key  = iflag_split(j,i)
<a name="l00120"></a>00120              call random_number(ran)      <span class="comment">! random set</span>
<a name="l00121"></a>00121              ran = 0.01d0 + 0.99d0*ran
<a name="l00122"></a>00122 <span class="comment">!            ran  = 0.5d0                 ! center of cartesian</span>
<a name="l00123"></a>00123              dv01 =  split(key-1) + ran*(split(key)-split(key-1))
<a name="l00124"></a>00124              sample(j,i) = bound(1,j) + (bound(2,j)-bound(1,j))*dv01
<a name="l00125"></a>00125 410       continue
<a name="l00126"></a>00126 400     continue
<a name="l00127"></a>00127 
<a name="l00128"></a>00128 <span class="keyword">        end subroutine Latin_Hypercube</span>
<a name="l00129"></a>00129 <span class="keyword">end module</span>
</pre></div></div>
</div>
<!-- window showing the filter options -->
<div id="MSearchSelectWindow"
     onmouseover="return searchBox.OnSearchSelectShow()"
     onmouseout="return searchBox.OnSearchSelectHide()"
     onkeydown="return searchBox.OnSearchSelectKey(event)">
<a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(0)"><span class="SelectionMark">&#160;</span>All</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(1)"><span class="SelectionMark">&#160;</span>Classes</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(2)"><span class="SelectionMark">&#160;</span>Files</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(3)"><span class="SelectionMark">&#160;</span>Functions</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(4)"><span class="SelectionMark">&#160;</span>Variables</a></div>

<!-- iframe showing the search results (closed by default) -->
<div id="MSearchResultsWindow">
<iframe src="javascript:void(0)" frameborder="0" 
        name="MSearchResults" id="MSearchResults">
</iframe>
</div>

<hr class="footer"/><address class="footer"><small>
Copyright (C) 2012 Brian A. Lockwood <br>
Generated on Tue May 1 2012 17:15:55 for Kriging Library
by&#160;<a href="http://www.doxygen.org/index.html"><img class="footer"
							 src="doxygen.png"
							 alt="doxygen"/></a>
1.7.3
</small></address>
